library(Seurat)
library(ggplot2)
library(cowplot)
library(scater)
library(scran)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(dplyr)
library(Matrix)
library(clustree)
library(pheatmap)
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")

#单组分析####
#**********************************************
# 需要把数据放到不同文件目录下，且文件名应该是
# barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz
# 没有gz后缀也可以
#**********************************************

##ADJ1
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ1")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ1.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
pbmc
saveRDS(pbmc, file = "./data/temp/ADJ1.rds")


##ADJ2
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ2")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ2.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/ADJ2.rds")

##ADJ3
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ3")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ3.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/ADJ3.rds")

##ADJ4
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ4")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ4.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/ADJ4.rds")

##ADJ5
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ5")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ5.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/ADJ5.rds")

##ADJ6
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ6")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/ADJ6.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/ADJ6.rds")

##PDAC1
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC1")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC1.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC1.rds")


##PDAC2
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC2")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC2.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC2.rds")

##PDAC3
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC3")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC3.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC3.rds")

##PDAC4
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC4")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC4.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC4.rds")


##PDAC5
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC5")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC5.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC5.rds")



##PDAC6
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/PDAC6")
pbmc <- CreateSeuratObject(counts = a1, project = "GSE212966", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="./data/output/PDAC6.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "./data/temp/PDAC6.rds")


# 总体分析####
##1. 读取分组数据####
# 使用for循环创建变量并读入数据
num_vars <- 6
for (i in 1:num_vars) {
  # if(i>=3 && i<=5){next}  #ADJ3，ADJ4，ADJ5数据缺失，跳过执行
  var_name <- paste("ADJ", i, sep = "")
  assign(var_name,readRDS(file=paste0("./data/temp/",var_name,".rds")))# 为变量赋值
}
for (i in 1:num_vars) {
  var_name <- paste("PDAC", i, sep = "")
  assign(var_name,readRDS(file=paste0("./data/temp/",var_name,".rds")))# 为变量赋值
}

##2. 标注分组信息####
#ADJ
ADJ1<-RenameCells(ADJ1,add.cell.id="ADJ1",for.merge=T)
ADJ1@meta.data$tech<-"Normal"
ADJ1@meta.data$celltype<-"Normal_ADJ1"

ADJ2<-RenameCells(ADJ2,add.cell.id="ADJ2",for.merge=T)
ADJ2@meta.data$tech<-"Normal"
ADJ2@meta.data$celltype<-"Normal_ADJ2"

ADJ3<-RenameCells(ADJ3,add.cell.id="ADJ3",for.merge=T)
ADJ3@meta.data$tech<-"Normal"
ADJ3@meta.data$celltype<-"Normal_ADJ3"

ADJ4<-RenameCells(ADJ4,add.cell.id="ADJ4",for.merge=T)
ADJ4@meta.data$tech<-"Normal"
ADJ4@meta.data$celltype<-"Normal_ADJ4"

ADJ5<-RenameCells(ADJ5,add.cell.id="ADJ5",for.merge=T)
ADJ5@meta.data$tech<-"Normal"
ADJ5@meta.data$celltype<-"Normal_ADJ5"

ADJ6<-RenameCells(ADJ6,add.cell.id="ADJ6",for.merge=T)
ADJ6@meta.data$tech<-"Normal"
ADJ6@meta.data$celltype<-"Normal_ADJ6"

#PDAC
PDAC1<-RenameCells(PDAC1,add.cell.id="PDAC1",for.merge=T)
PDAC1@meta.data$tech<-"Tumor"
PDAC1@meta.data$celltype<-"Tumor_PDAC1"

PDAC2<-RenameCells(PDAC2,add.cell.id="PDAC2",for.merge=T)
PDAC2@meta.data$tech<-"Tumor"
PDAC2@meta.data$celltype<-"Tumor_PDAC2"

PDAC3<-RenameCells(PDAC3,add.cell.id="PDAC3",for.merge=T)
PDAC3@meta.data$tech<-"Tumor"
PDAC3@meta.data$celltype<-"Tumor_PDAC3"

PDAC4<-RenameCells(PDAC4,add.cell.id="PDAC4",for.merge=T)
PDAC4@meta.data$tech<-"Tumor"
PDAC4@meta.data$celltype<-"TUmor_PDAC4"

PDAC5<-RenameCells(PDAC5,add.cell.id="PDAC5",for.merge=T)
PDAC5@meta.data$tech<-"Tumor"
PDAC5@meta.data$celltype<-"Tumor_PDAC5"

PDAC6<-RenameCells(PDAC6,add.cell.id="PDAC6",for.merge=T)
PDAC6@meta.data$tech<-"Tumor"
PDAC6@meta.data$celltype<-"Tumor_PDAC6"


##3. 数据合并####
#这里要把Normal放在后面
ALL_ADJ <- merge(ADJ1, y = c(ADJ2, ADJ3, ADJ4, ADJ5, ADJ6))
ALL_PDAC <- merge(PDAC1, y = c(PDAC2, PDAC3, PDAC4, PDAC5, PDAC6))
ALL <- merge(ALL_PDAC, ALL_ADJ)
saveRDS(ALL,file="./data/temp/ALL_before_integrate.rds")
# before integrate
library(Seurat)
library(ggplot2)

hms <- ALL
hms<-readRDS(file="./data/temp/ALL_before_integrate.rds")

hms[["percent.mt"]] <- PercentageFeatureSet(hms, pattern = "^MT-")
dim(hms[["percent.mt"]])#全为0的话是不对，修改上面的pattern正则表达
p1 <- VlnPlot(hms, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0)
p1
p2 <- FeatureScatter(hms, feature1 = "nCount_RNA", feature2 = "percent.mt")
p2
p3 <- FeatureScatter(hms, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p3
pdf(file="./data/output/ALL_Feature_Count_Percent_VlinPlot.pdf", width = 13, height = 6)
p1
p2 + p3
dev.off()

# 标准化+降维
pancreas <- NormalizeData(object = hms, normalization.method = "LogNormalize", scale.factor = 1e4)
pancreas <- FindVariableFeatures(pancreas, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
pancreas <- ScaleData(pancreas, verbose = FALSE)
pancreas <- RunPCA(pancreas, npcs = 30, verbose = FALSE)
pancreas <- RunUMAP(pancreas, reduction = "pca", dims = 1:30)
p1 <- DimPlot(pancreas, reduction = "umap", group.by = "tech")
p1
p2 <- DimPlot(pancreas, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE)
p2
pdf(file="./data/output/ALL_before_integrate_DimPlot.pdf", width = 15, height = 6)
plot_grid(p1,p2)
dev.off()

# after integrate
pancreas.list <- SplitObject(pancreas, split.by = "celltype")
for (i in 1: length(pancreas.list)) {
  pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
  pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
####🪲🐞🪲🐞🪲🐞####
#这里务必更改顺序，把参考量基线量放在后面，差异分析的时候比的内容就反
reference.list <- pancreas.list[c("Tumor_PDAC1","Tumor_PDAC2", "Tumor_PDAC3", "Tumor_PDAC4", "Tumor_PDAC5", "Tumor_PDAC6", 
                                  "Normal_ADJ1","Normal_ADJ2","Normal_ADJ3","Normal_ADJ4","Normal_ADJ5","Normal_ADJ6",)]

#🌟FindIntegrationAnchors函数费时费内存，运行前保存好数据，释放内存
saveRDS(reference.list,file="./data/temp/reference_list.rds")
rm(list=setdiff(ls(), "reference.list"))#释放其余内存
reference.list <- readRDS(file="./data/temp/reference_list.rds")

##4. 整合分析####
pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

#🌟IntegrateData函数费时费内存，运行前保存好数据，释放内存
saveRDS(pancreas.anchors,file="./data/temp/pancreas_anchors.rds")
rm(list=setdiff(ls(), "pancreas.anchors"))#释放其余内存
pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)
saveRDS(pancreas.integrated,file="./data/temp/pancreas_integrated1.rds")

#设置默认Assay
DefaultAssay(pancreas.integrated) <- "integrated"
#ScaleData对数据缩放和居中
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)
###PCA降维####
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)
print(pancreas.integrated[["pca"]], dims = 1:30, nfeatures = 5)
pdf(file="./data/output/ALL_PCA降维_after_integrated_DimPlot.pdf", width = 13, height = 6)
p1 <- VizDimLoadings(pancreas.integrated, dims = 1:2, reduction = "pca")
p2 <- DimPlot(pancreas.integrated, reduction = "pca")
plot_grid(p1,p2)
dev.off()
###UMAP降维####
pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30)
pdf(file="./data/output/ALL_UMAP降维_after_integrated_DimPlot.pdf", width = 13, height = 6)
p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype")
plot_grid(p1,p2)
dev.off()
saveRDS(pancreas.integrated, file = "./data/temp/hms_after_integrated.rds")
# saveRDS(pancreas.integrated, file = "./data/temp/DAC_PDAC_after_integrated.rds")

#聚类获得亚群####

hms_individual_integrated <- pancreas.integrated
hms_individual_integrated<-readRDS(file="./data/temp/hms_after_integrated.rds")
# p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype")
# p1

##1. 确定参数dim⚠️####
#费时，可以不跑
pbmc <- JackStraw(hms_individual_integrated, num.replicate = 100,dims = 40)
pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
pdf("./data/output/ALL_JackStrawPlot_ElbowPlot.pdf",width = 13,height = 6)
JackStrawPlot(pbmc, dims = 1:30)
ElbowPlot(pbmc,ndims = 30)
# find how many 15 cluster
ElbowPlot(hms_individual_integrated)
dev.off()

##2. 亚群聚类####
###2.1 亚组分群####
####Clustree预测分群####
hms_neighbor<- FindNeighbors(hms_individual_integrated, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
pdf("./data/output/ALL_Clustree.pdf",width = 13,height = 12)
clustree(obj)
dev.off()

#开始分群
rm(list=ls())
hms_individual_integrated<-readRDS(file="./data/temp/hms_after_integrated.rds")
hms_neighbor<- FindNeighbors(hms_individual_integrated, dims = 1:30)

####resolution = 1.2####
#1.2为推荐最大值
hms_cluster <- FindClusters(hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
saveRDS(hms_cluster, file = "./data/temp/hms_cluster_test_1.2.rds")
#未标注UMAP图
pdf("./data/output/ALL_DimPlot_FindClusters_1.2.pdf",width = 13,height = 12)
DimPlot(hms_cluster, reduction = "umap", label = TRUE)# pt.size = 0.5
dev.off()

####resolution = 4⚠️####
#，测试值，对比效果（不执行）
#分出了77个亚群，感觉不太必要大分类这么细，亚分类在细分时候应该是可以的
hms_cluster_4 <- FindClusters(hms_neighbor, resolution = 4)
head(Idents(hms_cluster_4), 5)
hms_cluster_4<- RunUMAP(hms_cluster_4, dims = 1:30)

pdf("./data/output/ALL_DimPlot_FindClusters_4.pdf",width = 13,height = 12)
DimPlot(hms_cluster_4, reduction = "umap")
dev.off()
saveRDS(hms_cluster, file = "./data/temp/hms_cluster_test_4.rds")


###2.2 确定markers####
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/output/ALL_cellMarkers_1.2.txt",sep="\t",row.names=F,quote=F)

####top20_cell_markers####
#挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="./data/output/ALL_top20_cell_markers_1.2.csv",top20)
p1 <- DoHeatmap(hms_cluster, features = top20$gene) + NoLegend()
pdf("./data/output/ALL_top20_cell_DoHeatmap.pdf",width = 13,height = 12)
p1
dev.off()
#整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/output/ALL_top20_marker_genes_1.2.csv",top20_table,row.names=F)

#查看信息
slotNames(hms_cluster) # 细胞及细胞中基因与RNA数量
hms_cluster@assays     # assay
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)

hms_cluster<-readRDS(file="./data/temp/hms_cluster_test_1.2.rds")
DimPlot(hms_cluster, reduction = "umap")

##3. SingleR自动注释(⚠️仅做参考)####
# SingleR https://zhuanlan.zhihu.com/p/448269413
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("SingleR")
library(SingleR)
library(Seurat)
library(pheatmap)
# load("HumanPrimaryCellAtlas_hpca.se_human.RData")
load("./data/input/BlueprintEncode_bpe.se_human.RData")

#进行SingleR注释
hms_cluster_for_SingleR <- GetAssayData(hms_cluster, slot="data") ##获取标准化矩阵
hms_cluster.hesc <- SingleR(test = hms_cluster_for_SingleR, ref = bpe.se, labels = bpe.se$label.main) #
hms_cluster.hesc
hms_cluster.hesc$labels

#seurat 和 SingleR的table表
a <- table(hms_cluster.hesc$labels, hms_cluster@meta.data$seurat_clusters)# pbmc@meta.data  pbmc的meta文件，包含了seurat的聚类结果
write.csv(a, "./data/output/SingleR_by_pbe.csv")

#保存亚群名到hms_cluster
hms_cluster@meta.data$labels <-hms_cluster.hesc$labels
saveRDS(hms_cluster, file = "./data/temp/hms_cluster_by_SingleR_labels.rds")

pdf("./data/output/ALL_Dimplot_SingleR_by_bpe.pdf",width = 21,height = 8)
DimPlot(hms_cluster, group.by = c("seurat_clusters", "labels"),reduction = "umap",label = T)
dev.off()

##4. 注释基因####
hms_cluster<-readRDS(file="./data/temp/hms_cluster_test_1.2.rds")
#⚠️第二次的标注分群
# new.cluster.ids <- c("T_cell"	,"T_cell"	,"B_cell"	,"Fibroblast"	,"T_cell"	,#4
#                      "Fibroblast"	,"Macrophage"	,"Endothelial_cell"	,"Ductal_cell"	,"T_cell"	,#9
#                      "Schwann_cell"	,"Epithelial_cell"	,"Neutrophil"	,"Epithelial_cell"	,"Epithelial_cell"	,#14
#                      "Stellate_cell"	,"Epithelial_cell"	,"T_cell"	,"Macrophage"	,"Fibroblast"	,#19
#                      "T_cell"	,"Acinar_cell"	,"Dendritic_cell"	,"Mast_cell"	,"T_cell"	,#24
#                      "Acinar_cell"	,"Plasma_cell"	,"T_cell"	,"Epithelial_cell"	,"Plasma_cell"	,#29
#                      "Endocrine_cell"	,"T_cell"	,"Schwann_cell"	,"Endothelial_cell"	,"Epithelial_cell"	,#34
#                      "Epithelial_cell"	,"T_cell"
# )
#⭐20240513 第三次标注分群 将epithelial 改为ductal
new.cluster.ids <- c("T_cell"	,"T_cell"	,"B_cell"	,"Fibroblast"	,"T_cell"	,#4
                     "Fibroblast"	,"Macrophage"	,"Endothelial_cell"	,"Duct_epithelial_cell"	,"T_cell"	,#9
                     "Schwann_cell"	,"Duct_epithelial_cell"	,"Neutrophil"	,"Duct_epithelial_cell"	,"Duct_epithelial_cell"	,#14
                     "Stellate_cell"	,"Duct_epithelial_cell"	,"T_cell"	,"Macrophage"	,"Fibroblast"	,#19
                     "T_cell"	,"Acinar_cell"	,"Dendritic_cell"	,"Mast_cell"	,"T_cell"	,#24
                     "Acinar_cell"	,"Plasma_cell"	,"T_cell"	,"Duct_epithelial_cell"	,"Plasma_cell"	,#29
                     "Endocrine_cell"	,"T_cell"	,"Schwann_cell"	,"Endothelial_cell"	,"Duct_epithelial_cell"	,#34
                     "Duct_epithelial_cell"	,"T_cell"
)
levels(hms_cluster)
names(new.cluster.ids) <- levels(hms_cluster)
#在hms_cluster@active.ident添加分类标记
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
hms_cluster_id@active.ident
# hms_cluster_id <- hms_cluster
# 调整active.idents顺序，方便图注顺序
levels(hms_cluster_id)
levels(hms_cluster_id) <- c("Acinar_cell", "B_cell",  "Dendritic_cell"  ,"Duct_epithelial_cell", "Endocrine_cell", 
                   "Endothelial_cell", "Fibroblast", "Macrophage", "Mast_cell", 
                   "Neutrophil", "Plasma_cell", "Schwann_cell", "Stellate_cell","T_cell")
levels(hms_cluster_id)
# save
saveRDS(hms_cluster_id, file = "./data/temp/hms_cluster_id_test1.3.rds")

### UMAP出图####
pdf("./data/output/ALL_DimPlot_已标注.pdf",width =22 ,height = 10)
dp1 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5 )
dp2 <- DimPlot(hms_cluster_id, reduction = "umap", label = FALSE, pt.size = 0.5 )
dp1+dp2
dev.off()

# 按肿瘤组和正常组分别展示
pdf("./data/output/ALL_DimPlot_已标注_Normal_Tumor.pdf",width = 19,height = 10)
DimPlot(hms_cluster_id, reduction = "umap", split.by = "tech",label = TRUE, pt.size = 0.5 )
dev.off()


##5. 注释验证####
rm(list=setdiff(ls(), "hms_cluster_id"))#释放其余内存
# OR
hms_cluster_id <- readRDS("./data/temp/hms_cluster_id_test.rds")

###5.1 方法1: DotPlot####
markers1 <- c(
  "PRSS1","PRSS2","CTRB1",# Acinar cell"REG1B",
  "MS4A1", "CD79A","LTB",# B cell "CD19", "CD83",
  "EREG", "IL1B", "HLA-DPA1", # Dendritic_cell
  "KRT7", "KRT17", "KRT19",  # Ductal cell"MMP7",
  "CHGA", "CHGB", "INS",  # Endocrine cell"TTR", 
  "VWF",  "PECAM1", "CLDN5",# Endothelial cell"KDR",
  #"KRT17", "S100A6", "TFF1",#Epithelial_cell
  "COL1A1", "COL1A2", "COL3A1",# Fibroblast "CTHRC1",
  "C1QA", "C1QB" ,"FTL",  # Macrophage
  "TPSAB1","TPSB2", "CPA3", #Mast_cell"KIT", 
  #"GNLY", "GZMB", "KLRD1",# NK cells "PRF1",
  "G0S2", "S100A9", "S100A8",# Neitrophil "CXCL8",
  "IGHA1", "JCHAIN", "MZB1",# Plasma cell
  "CDH19", "ITGB8","GPM6B",# Schwann cell,"SCN7A"
  "ADIRF", "RGS5", "PDGFRB",  # Stellate cell "COL4A1", 
  "CCL5","GNLY", "GZMK")   # T cell 待定"BATF", 
dp1 <-DotPlot(hms_cluster_id, features = markers1, dot.scale = 5) + RotatedAxis()
dp1
#####DotPlot_全聚类_标注验证####
pdf("./data/output/All_DotPlot_标注验证.pdf",width = 13,height = 6)
dp1
dev.off()


###5.2 方法2: FeaturePlot####
fp1 <- FeaturePlot(hms_cluster_id,reduction ="umap",features =c("HBB"))
fp1
###5.3 方法3: VlnPlot####
fp2 <- VlnPlot(hms_cluster_id, features = c("MS4A1", "CD79A"))
fp2

##6. 分群绘图####
hms_cluster_id<-readRDS("./data/temp/hms_cluster_id_test.rds")

Acinar<-subset(hms_cluster_id, idents=c('Acinar_cell'))
pdf("./data/output/Acinar_DimPlot.pdf",width = 8,height = 6)
DimPlot(Acinar, reduction = "umap")
dev.off()
saveRDS(Acinar, file="./data/temp/Acinar.rds")

B<-subset(hms_cluster_id, idents=c('B_cell'))
pdf("./data/output/B_DimPlot.pdf",width = 8,height = 6)
DimPlot(B, reduction = "umap")
dev.off()
saveRDS(B, file="./data/temp/B.rds")

Dendritic<-subset(hms_cluster_id, idents=c('Dendritic_cell'))
pdf("./data/output/Dendritic_DimPlot.pdf",width = 8,height = 6)
DimPlot(Dendritic, reduction = "umap")
dev.off()
saveRDS(Dendritic, file="./data/temp/Dendritic.rds")

Ductal<-subset(hms_cluster_id, idents=c('Ductal_cell'))
pdf("./data/output/Ductal_DimPlot.pdf",width = 8,height = 6)
DimPlot(Ductal, reduction = "umap")
dev.off()
saveRDS(Ductal, file="./data/temp/Ductal.rds")

Endocrine<-subset(hms_cluster_id, idents=c('Endocrine_cell'))
pdf("./data/output/Endocrine_DimPlot.pdf",width = 8,height = 6)
DimPlot(Endocrine, reduction = "umap")
dev.off()
saveRDS(Endocrine, file="./data/temp/Endocrine.rds")

Endothelial<-subset(hms_cluster_id, idents=c('Endothelial_cell'))
pdf("./data/output/Endothelial_DimPlot.pdf",width = 8,height = 6)
DimPlot(Endothelial, reduction = "umap")
dev.off()
saveRDS(Endothelial, file="./data/temp/Endothelial.rds")

Epithelial<-subset(hms_cluster_id, idents=c('Epithelial_cell'))
pdf("./data/output/Epithelial_DimPlot.pdf",width = 8,height = 6)
DimPlot(Epithelial, reduction = "umap")
dev.off()
saveRDS(Epithelial, file="./data/temp/Epithelial.rds")

Fibroblast<-subset(hms_cluster_id, idents=c('Fibroblast'))
pdf("./data/output/Fibroblast_DimPlot.pdf",width = 8,height = 6)
DimPlot(Fibroblast, reduction = "umap")
dev.off()
saveRDS(Fibroblast, file="./data/temp/Fibroblast.rds")

Macrophage<-subset(hms_cluster_id, idents=c('Macrophage'))
pdf("./data/output/Macrophage_DimPlot.pdf",width = 8,height = 6)
DimPlot(Macrophage, reduction = "umap")
dev.off()
saveRDS(Macrophage, file="./data/temp/Macrophage.rds")

Mast<-subset(hms_cluster_id, idents=c('Mast_cell'))
pdf("./data/output/Mast_DimPlot.pdf",width = 8,height = 6)
DimPlot(Mast, reduction = "umap")
dev.off()
saveRDS(Mast, file="./data/temp/Mast.rds")

Neutrophil<-subset(hms_cluster_id, idents=c('Neutrophil'))
pdf("./data/output/Neutrophil_DimPlot.pdf",width = 8,height = 6)
DimPlot(Neutrophil, reduction = "umap")
dev.off()
saveRDS(Neutrophil, file="./data/temp/Neutrophil.rds")

Mast<-subset(hms_cluster_id, idents=c('Mast_cell'))
pdf("./data/output/Mast_DimPlot.pdf",width = 8,height = 6)
DimPlot(Mast, reduction = "umap")
dev.off()
saveRDS(Mast, file="./data/temp/Mast.rds")

Neutrophil<-subset(hms_cluster_id, idents=c('Neutrophil'))
pdf("./data/output/Neutrophil_DimPlot.pdf",width = 8,height = 6)
DimPlot(Neutrophil, reduction = "umap")
dev.off()
saveRDS(Neutrophil, file="./data/temp/Neutrophil.rds")

Plasma<-subset(hms_cluster_id, idents=c('Plasma_cell'))
pdf("./data/output/Plasma_DimPlot.pdf",width = 8,height = 6)
DimPlot(Plasma, reduction = "umap")
dev.off()
saveRDS(Plasma, file="./data/temp/Plasma.rds")

Schwann<-subset(hms_cluster_id, idents=c('Schwann_cell'))
pdf("./data/output/Schwann_DimPlot.pdf",width = 8,height = 6)
DimPlot(Schwann, reduction = "umap")
dev.off()
saveRDS(Schwann, file="./data/temp/Schwann.rds")

Stellate<-subset(hms_cluster_id, idents=c('Stellate_cell'))
pdf("./data/output/Stellate_DimPlot.pdf",width = 8,height = 6)
DimPlot(Stellate, reduction = "umap")
dev.off()
saveRDS(Stellate, file="./data/temp/Stellate.rds")

T<-subset(hms_cluster_id, idents=c('T_cell'))
pdf("./data/output/T_DimPlot.pdf",width = 8,height = 6)
DimPlot(T, reduction = "umap")
dev.off()
saveRDS(T, file="./data/temp/T.rds")

###展示亚群对比####
##input each cluster
Acinar <- readRDS(file = "./data/temp/Acinar.rds")
B <- readRDS(file = "./data/temp/B.rds")
Dendritic <- readRDS(file = "./data/temp/Dendritic.rds")
Ductal <- readRDS(file = "./data/temp/Ductal.rds")
Endocrine <- readRDS(file = "./data/temp/Endocrine.rds")
Endothelial <- readRDS(file = "./data/temp/Endothelial.rds")
Epithelial <- readRDS(file = "./data/temp/Epithelial.rds")
Fibroblast <- readRDS(file = "./data/temp/Fibroblast.rds")
Macrophage <- readRDS(file = "./data/temp/Macrophage.rds")
# NK <- readRDS(file = "./data/temp/NK.rds")
Mast <- readRDS(file = "./data/temp/Mast.rds")
Neutrophil <- readRDS(file = "./data/temp/NK.rds")
Plasma <- readRDS(file = "./data/temp/Plasma.rds")
Schwann <- readRDS(file = "./data/temp/Schwann.rds")
Stellate <- readRDS(file = "./data/temp/Stellate.rds")
T <- readRDS(file = "./data/temp/T.rds")

pdf("./data/output/SUB_DimPlot_Normal_Tumor合并.pdf",width = 13,height = 6)
DimPlot(Acinar, reduction = "umap", split.by = "tech")
DimPlot(B, reduction = "umap", split.by = "tech")
DimPlot(Dendritic, reduction = "umap", split.by = "tech")
DimPlot(Ductal, reduction = "umap", split.by = "tech")
DimPlot(Endocrine, reduction = "umap", split.by = "tech")
DimPlot(Endothelial, reduction = "umap", split.by = "tech")
DimPlot(Epithelial, reduction = "umap", split.by = "tech")
DimPlot(Fibroblast, reduction = "umap", split.by = "tech")
DimPlot(Macrophage, reduction = "umap", split.by = "tech")
DimPlot(Mast, reduction = "umap", split.by = "tech")
# DimPlot(NK, reduction = "umap", split.by = "tech")
DimPlot(Neutrophil, reduction = "umap", split.by = "tech")
DimPlot(Plasma, reduction = "umap", split.by = "tech")
DimPlot(Schwann, reduction = "umap", split.by = "tech")
DimPlot(Stellate, reduction = "umap", split.by = "tech")
DimPlot(T, reduction = "umap", split.by = "tech")
dev.off()


#亚群再聚类####
## 1. T亚群####
###1.1 亚群聚类####
time_start <- Sys.time()
T <- readRDS("./data/temp/T.rds")
# # 可以不跑 JackStraw，ScoreJackStraw
# pbmc <- JackStraw(T, num.replicate = 100,dims = 30)#40
# pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
# pdf("./data/output/T_JackStrawPlot_ElbowPlot.pdf",width = 8,height = 6)
# JackStrawPlot(pbmc, dims = 1:30)
# ElbowPlot(pbmc,ndims = 30)
# dev.off()
# # find how many 15 cluster
# hms_neighbor<- FindNeighbors(pbmc, dims = 1:30)
hms_neighbor<- FindNeighbors(T, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
pdf("./data/output/T_clustree.pdf",width = 8,height = 6)
clustree(obj)
dev.off()

####resolution = 1.2####
hms_cluster <- FindClusters(hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
#UMAP可视化
pdf("./data/output/T_未标注_DimPlot.pdf",width = 7,height = 6)
DimPlot(hms_cluster, reduction = "umap",label = TRUE, pt.size = 0.5)
dev.off()
# save
saveRDS(hms_cluster, file = "./data/temp/T_test_1.2.rds")

# cellmarkers用于后续标注
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/temp/T_cellMarkers.txt",sep="\t",row.names=F,quote=F)
# 挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="./data/temp/T_cellmarkers_top20.csv",top20)
# 整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/T_cellmarkers_top20_genes_1.2.csv",top20_table,row.names=F)

#耗时
time_end <- Sys.time()
time_cost <- time_end-time_start
time_cost
# Time difference of 9.383156 mins

####resolution = 4####
time_start <- Sys.time()
T <- readRDS("./data/temp/T.rds")
hms_neighbor_4<- FindNeighbors(T, dims = 1:30)
hms_cluster_4 <- FindClusters(hms_neighbor_4, resolution = 4)
head(Idents(hms_cluster_4), 5)
hms_cluster_4<- RunUMAP(hms_cluster_4, dims = 1:30)
#UMAP可视化
pdf("./data/output/T_未标注_DimPlot_4.pdf",width = 7,height = 6)
DimPlot(hms_cluster_4, reduction = "umap",label = TRUE, pt.size = 0.5)
dev.off()
# save
saveRDS(hms_cluster_4, file = "./data/temp/T_test_4.rds")

# cellmarkers用于后续标注
scRNA.markers_4 <- FindAllMarkers(hms_cluster_4, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers_4,file="./data/temp/T_cellMarkers_4.txt",sep="\t",row.names=F,quote=F)
# 挑选每个细胞亚群中特意高表达的20个基因
top20_4 <- scRNA.markers_4 %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="./data/temp/T_cellmarkers_top20_4.csv",top20)
# 整理成表格，只显示基因名字
top20_table_4=unstack(top20_4, gene ~ cluster)
names(top20_table_4)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/T_cellmarkers_top20_genes_4.csv",top20_table_4,row.names=F)

#耗时
time_end <- Sys.time()
time_cost <- time_end-time_start
time_cost
# Time difference of 12.90657 mins

###1.2 亚群注释####
hms_cluster<-readRDS("./data/temp/T_test_1.2.rds")
DimPlot(hms_cluster, reduction = "umap", label = TRUE)

#细胞及细胞中基因与RNA数量
slotNames(hms_cluster)
#assay
hms_cluster@assays
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)

new.cluster.ids <- c("CD4_Tn", "CD4_Th", "CD8_Te", "CD4_Tem", "CD8_Tc", 
                     "CD8_Tem", "Endocrine_cell", "NKT", "CD8_Trm", "CD8_Tex", 
                     "NKT", "CD4_Treg", "Endocrine_cell", "Acinar_cell", "CD4_Treg", 
                     "CD8_Tc", "Macrophage", "CD4_Treg", "Dendritic_cell", "CD4_Tem")
names(new.cluster.ids) <- levels(hms_cluster)
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
# levels()修改ident的顺序，解决后期作图排序
levels(hms_cluster_id)
levels(hms_cluster_id) <- c("CD4_Tem", "CD4_Th", "CD4_Tn", "CD4_Treg", 
                            "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Tn", "CD8_Trm", 
                            "Acinar_cell","Dendritic_cell", "Endocrine_cell", "Macrophage", "NKT")
levels(hms_cluster_id)
DimPlot(hms_cluster_id, reduction = "umap", label = TRUE)
# save
# saveRDS(hms_cluster_id, file = "./data/temp/T_cluster_id_test.rds")#2024.1.31,CD8分群修改
saveRDS(hms_cluster_id, file = "./data/temp/T_cluster_id_test_1.rds")


###1.3 注释验证####
#绘图准备： 读入+亚群提取
# hms_cluster_id<-readRDS("./data/temp/T_cluster_id_test.rds")
hms_cluster_id<-readRDS("./data/temp/T_cluster_id_test_1.rds")
sub_cluster_cell <- WhichCells(hms_cluster_id, idents = c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
                                                        "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem","CD8_Tex", "CD8_Trm"))
sub_cluster_id <- subset(hms_cluster_id,idents=c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
                                               "NKT","CD8_Tc","CD8_Te","CD8_Tem","CD8_Tex", "CD8_Trm"))

####1.3.1 UMAP_Dimplot：某个基因分布UMAP####
#单独的基因分布,可以选择不同基因显示，按照需要作图
#可能后面需要写个循环出图
p1 <- FeaturePlot(hms_cluster_id,reduction ="umap",features =c("GNLY"))
p1
####1.3.2 Marker_VlnPlot：某个基因在簇中的分布####
p2 <- VlnPlot(hms_cluster_id, features = c("GNLY"),pt.size = 0)
p2

pdf("./data/output/T_sub_标注验证1_1.pdf",width =10 ,height = 9.5)
p1
p2
dev.off()

#⭐未使用，单独设置颜色出图
# colors <- c("#942d8d","#3cb346","#00abf0","#d75427",
#             "#2e409a","#FB8072","#eeb401")#,"#7BAFDE"
# DimPlot(hms_cluster_id,label=T,cells = g1_treat, cols =colors )

# T亚群+其他
p3 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5) # label = FALSE
p3
# 仅T亚群
p4 <- DimPlot(hms_cluster_id, cells = sub_cluster_cell, reduction = "umap", label = TRUE,  pt.size = 0.5)
p4
pdf("./data/output/T_sub_标注验证2_1.pdf",width =10 ,height = 8.3)
p3
p4
dev.off()

# 仅T亚群对比：Normal_Tumor
p5 <- DimPlot(sub_cluster_id, cells = sub_cluster_cell,  reduction = "umap", split.by = "tech",label = TRUE, pt.size = 0.5 )
p5
p6 <- DimPlot(hms_cluster_id, reduction = "umap", group.by = "tech",repel =TRUE)
p6
pdf("./data/output/T_sub_标注验证3_1.pdf",width = 19,height = 10)
p5
dev.off()
pdf("./data/output/T_sub_标注验证4_1.pdf",width = 10,height = 9)
p6
dev.off()

####1.3.3 DotPlot: 全亚群-气泡热图####
#⚠️T亚群注释的并不清晰
#⚠️⚠️⚠️ene还要重新选的,可以再加gene
markers1 <- c(
  #"CD3D",
  "AQP3", "KLRB1", # CD4_Tem
  "CXCL13", "LTB",     # CD4_Th 
  "GPR183", "CCR7", "CD4",   # CD4_Tn
  "RGS1", "BATF", "TNFRSF4","CTLA4",     # CD4_Treg
  "DUSP2", "CCL5",       # CD8_Tc
   "GZMA",        # CD8_Te
  "GZMK", "CRTAM",       # CD8_Tem 
  "CCL4", "CCL4L2",      # CD8_Tex
  "CD69",  "GZMH",              # CD8_Trm 
  "GZMB", "GNLY"         # NKT
  )
p7 <-DotPlot(subset(hms_cluster_id,idents=c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
                                             "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")), #
              features = markers1, 
              dot.scale = 5) + RotatedAxis()
p7

#T亚群中非T细胞
markers2=c(
  "PRSS1","CTRB1","PRSS2", # Acinar cell
  "HLA-DRA", "HLA-DPA1","CD74",  # Dendritic cell
  "PCSK1N", "APOD", "ZG16B", # Endocrine_cell
  "TOP2A", "STMN1","MKI67") # Macrophages
p8 <-DotPlot(subset(hms_cluster_id,idents=c("Endocrine_cell", "Macrophage","Dendritic_cell", "Acinar_cell")), 
              features = markers2, 
              dot.scale = 5) + 
  RotatedAxis()
p8

pdf("./data/output/T_sub_标注验证5_气泡热图_1.pdf.pdf",width = 13,height = 6)
p7
p8
dev.off()

####1.3.4 RidgePlot：Cellmarker表达量-Normal Vs Tumo ####
#Cellmarker表达量:对比Normal Vs Tumo
p9 <- RidgePlot(sub_cluster_id, 
                features = 
                  c(
                    #"CD3D",
                    "AQP3", "KLRB1", # CD4_Tem
                    "CXCL13", "LTB",     # CD4_Th 
                    "GPR183", "CCR7", "CD4",   # CD4_Tn
                    "RGS1", "BATF", "TNFRSF4","CTLA4",     # CD4_Treg
                    "DUSP2", "CCL5",       # CD8_Tc
                    "GZMA",        # CD8_Te
                    "GZMK", "CRTAM",       # CD8_Tem 
                    "CCL4", "CCL4L2",      # CD8_Tex
                    "CD69",  "GZMH",              # CD8_Trm 
                    "GZMB", "GNLY"         # NKT
                  ),
                # c("AQP3", "KLRB1", # CD4_Tem
                #             "CXCL13","RGS1", # CD4_Th 
                #             "GPR183", "CCR7", # CD4_Tn
                #             "BATF", "TNFRSF4", # CD4_Treg
                #             "CCL5","DUSP2", # CD8_Tc
                #             "CCL4", "CCL4L2", # CD8_Te
                #             "CRTAM", "GZMK", # CD8_Tem 
                #             "GZMB", "GNLY",# NKT
                #             "GZMH", "GZMA") # CD8_Trm 
                cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())
p9
pdf("./data/output/T_sub_标注验证6_cellmarkers表达量对比_1.pdf",width = 10,height = 10)
p9
dev.off()

####1.3.5 DoHeatmap：全亚群-热图####
markers.to.plot <-
  # c("AQP3", "KLRB1", # CD4_Tem
  #                   "CXCL13","RGS1", # CD4_Th 
  #                   "GPR183", "CCR7", # CD4_Tn
  #                   "BATF", "TNFRSF4", # CD4_Treg
  #                   "CCL5","DUSP2", # CD8_Tc
  #                   "CCL4", "CCL4L2", # CD8_Te
  #                   "CRTAM", "GZMK", # CD8_Tem 
  #                   "GZMB", "GNLY",# NKT
  #                   "GZMH", "GZMA") # CD8_Trm 
c(
  #"CD3D",
  "AQP3", "KLRB1", # CD4_Tem
  "CXCL13", "LTB",     # CD4_Th 
  "GPR183", "CCR7", "CD4",   # CD4_Tn
  "RGS1", "BATF", "TNFRSF4","CTLA4",     # CD4_Treg
  "DUSP2", "CCL5",       # CD8_Tc
  "GZMA",        # CD8_Te
  "GZMK", "CRTAM",       # CD8_Tem 
  "CCL4", "CCL4L2",      # CD8_Tex
  "CD69",  "GZMH",              # CD8_Trm 
  "GZMB", "GNLY"         # NKT
)
p10 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
                features = markers.to.plot,
                size=5,
                group.bar.height=0.03)+
       scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p10
T_sub_heatmap<-p10$data
write.table(T_sub_heatmap,"./data/temp/T_sub_heatmap_1.csv")

p11 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
              features = markers.to.plot,
              size=5,
              group.by="celltype")+
       scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p11
T_sub_celltype_heatmap<-p11$data
write.table(T_sub_celltype_heatmap,"./data/temp/T_sub_celltype_heatmap_1.csv")

pdf("./data/output/T_sub_标注验证7_Heatmap_genes_celltype_1.pdf",width = 16,height = 13)
p10
p11
dev.off()







































































































































































































































##2. Fibroblast亚群####
Fibroblast <- readRDS("./data/temp/Fibroblast.rds")
# # 可以不跑 JackStraw，ScoreJackStraw
# pbmc <- JackStraw(T, num.replicate = 100,dims = 30)#40
# pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
# pdf("./data/output/T_JackStrawPlot_ElbowPlot.pdf",width = 8,height = 6)
# JackStrawPlot(pbmc, dims = 1:30)
# ElbowPlot(pbmc,ndims = 30)
# dev.off()
# # find how many 15 cluster
# hms_neighbor<- FindNeighbors(pbmc, dims = 1:30)
hms_neighbor<- FindNeighbors(Fibroblast, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
pdf("./data/output/Fibroblast_clustree.pdf",width = 8,height = 6)
clustree(obj)
dev.off()

####resolution = 1.2####
hms_cluster <- FindClusters(hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
#UMAP可视化
pdf("./data/output/Fibroblast_未标注_DimPlot.pdf",width = 7,height = 6)
DimPlot(hms_cluster, reduction = "umap",label = TRUE, pt.size = 0.5)
dev.off()
# save
saveRDS(hms_cluster, file = "./data/temp/Fibroblast_test_1.2.rds")

# cellmarkers用于后续标注
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/temp/Fibroblast_cellMarkers.txt",sep="\t",row.names=F,quote=F)
# 挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="./data/temp/Fibroblast_cellmarkers_top20.csv",top20)
# 整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/Fibroblast_cellmarkers_top20_genes_1.2.csv",top20_table,row.names=F)


 ###2.2 亚群注释####
hms_cluster<-readRDS("./data/temp2/Fibroblast_test_1.2.rds")
DimPlot(hms_cluster, reduction = "umap", label = TRUE)

#细胞及细胞中基因与RNA数量
slotNames(hms_cluster)
#assay
hms_cluster@assays
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)

new.cluster.ids <- c("myCAF"	,"iCAF"	,"iCAF"	,"iCAF"	,"iCAF"	,"myCAF"	,"mCAF"	,"iCAF"	,"iCAF"	,"myCAF"	,"myCAF"	,"iCAF"	,"myCAF"	,"myCAF"	,"myCAF"	,"mCAF"	,"apCAF")
names(new.cluster.ids) <- levels(hms_cluster)
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)


# levels()修改ident的顺序，解决后期作图排序
levels(hms_cluster_id)
levels(hms_cluster_id) <- c("apCAF", "iCAF", "mCAF", "myCAF","Epithelial_cell")
levels(hms_cluster_id)
# save
saveRDS(hms_cluster_id, file = "./data/temp/Fibroblast_cluster_id_test.rds")

###2.3 注释验证####
#绘图准备： 读入+亚群提取
hms_cluster_id<-readRDS("./data/temp/Fibroblast_cluster_id_test.rds")
sub_cluster_cell <- WhichCells(hms_cluster_id, idents = c("iCAF", "mCAF", "myCAF"))#"apCAF", 
sub_cluster_id <- subset(hms_cluster_id,idents=c("iCAF", "mCAF", "myCAF"))#"apCAF", 
#20240511 只保留三种亚群的数据⭐
levels(sub_cluster_id)
saveRDS(sub_cluster_id, file = "./data/temp/Fibroblast_3亚群_i-m-myCAF.rds")

####1.3.1 UMAP_Dimplot：某个基因分布UMAP####
#单独的基因分布,可以选择不同基因显示，按照需要作图
#可能后面需要写个循环出图
p1 <- FeaturePlot(hms_cluster_id,reduction ="umap",features =c("CD74"))
p1
####1.3.2 Marker_VlnPlot：某个基因在簇中的分布####
p2 <- VlnPlot(hms_cluster_id, features = c("CD74"),pt.size = 0)
p2

pdf("./data/output/Fibroblast_sub_标注验证1.pdf",width =10 ,height = 9.5)
p1
p2
dev.off()

#⭐未使用，单独设置颜色出图
# colors <- c("#942d8d","#3cb346","#00abf0","#d75427",
#             "#2e409a","#FB8072","#eeb401")#,"#7BAFDE"
# DimPlot(hms_cluster_id,label=T,cells = g1_treat, cols =colors )

# 亚群+其他
p3 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5) # label = FALSE
p3
# 仅T亚群
p4 <- DimPlot(hms_cluster_id, cells = sub_cluster_cell, reduction = "umap", label = TRUE,  pt.size = 0.5)
p4
pdf("./data/output/Fibroblast_sub_标注验证2.pdf",width =10 ,height = 8.3)
p3
p4
dev.off()

# 仅T亚群对比：Normal_Tumor
p5 <- DimPlot(sub_cluster_id, cells = sub_cluster_cell,  reduction = "umap",label = TRUE, pt.size = 0.5 )
p5
p6 <- DimPlot(sub_cluster_id, reduction = "umap",  split.by = "tech",repel =TRUE)
p6
pdf("./data/output/Fibroblast_分群UMAP_合并.pdf",width = 6,height = 5)
p5
dev.off()
pdf("./data/output/Fibroblast_分群UMAP_NT对比.pdf",width = 10,height = 5)
p6
dev.off()

####1.3.3 DotPlot: 全亚群-气泡热图####
#⚠️T亚群注释的并不清晰
#⚠️⚠️⚠️ene还要重新选的,可以再加gene
markers1 <- c(
  #"HLA-DRA","CD74","HLA-DRB1",# apCAF
  "FBLN1", "PLA2G2A","HSPA1A", # iCAF 
  "FN1", "COL1A2", "COL1A1", # mCAF
  "CST1", "CST2", "MYH11") # myCAF 
p7 <-DotPlot(subset(hms_cluster_id,idents= c( "iCAF", "mCAF", "myCAF")), #"apCAF",
             features = markers1, 
             dot.scale = 5) + RotatedAxis()
p7

#T亚群中非T细胞
markers2=c("FGFR2","CYP1B1","GCG")# Epithelial cell
p8 <-DotPlot(subset(hms_cluster_id,idents=c("Epithelial_cell")), 
             features = markers2, 
             dot.scale = 5) + 
  RotatedAxis()
p8

pdf("./data/output/Fibroblast_marker_气泡热图.pdf.pdf",width = 13,height = 6)
p7
p8
dev.off()

####1.3.4 RidgePlot：Cellmarker表达量-Normal Vs Tumo ####
#Cellmarker表达量:对比Normal Vs Tumo
p9 <- RidgePlot(sub_cluster_id, 
                features = c(
                  "HLA-DRA","CD74","HLA-DRB1",# apCAF
                  "FBLN1", "PLA2G2A","HSPA1A", # iCAF 
                  "FN1", "COL1A2", "COL1A1", # mCAF
                  "CST1", "CST2", "MYH11") , # myCAF 
                cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())
p9
pdf("./data/output/Fibroblast_sub_标注验证6_cellmarkers表达量对比.pdf",width = 10,height = 10)
p9
dev.off()

####1.3.5 DoHeatmap：全亚群-热图####
markers.to.plot <-c(
  "HLA-DRA","CD74","HLA-DRB1",# apCAF
  "FBLN1", "PLA2G2A","HSPA1A", # iCAF 
  "FN1", "COL1A2", "COL1A1", # mCAF
  "CST1", "CST2", "IGFBP2") # myCAF 
p10 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
                 features = markers.to.plot,
                 size=5,
                 group.bar.height=0.03)+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p10
T_sub_heatmap<-p10$data
write.table(T_sub_heatmap,"./data/temp/Fibroblast_sub_heatmap.csv")

p11 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
                 features = markers.to.plot,
                 size=5,
                 group.by="celltype")+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p11
T_sub_celltype_heatmap<-p11$data
write.table(T_sub_celltype_heatmap,"./data/temp/Fibroblast_sub_celltype_heatmap.csv")

pdf("./data/output/Fibroblast_sub_标注验证7_Heatmap_genes_celltype.pdf",width = 16,height = 13)
p10
p11
dev.off()

##3. Duct_epithelial亚群####
hms_cluster_id<-readRDS("./data/temp/hms_cluster_id_test1.3.rds")
Duct_epithelial<-subset(hms_cluster_id, idents=c('Duct_epithelial_cell'))
pdf("./data/output/Duct_epithelial_DimPlot_初分Normal_Tumor.pdf",width = 8,height = 6)
DimPlot(Duct_epithelial, reduction = "umap",split.by = "tech")
dev.off()
saveRDS(Duct_epithelial, file="./data/temp/Duct_epithelial.rds")

Duct_epithelial <- readRDS("./data/temp/Duct_epithelial.rds")
# # 可以不跑 JackStraw，ScoreJackStraw
# pbmc <- JackStraw(T, num.replicate = 100,dims = 30)#40
# pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
# pdf("./data/output/T_JackStrawPlot_ElbowPlot.pdf",width = 8,height = 6)
# JackStrawPlot(pbmc, dims = 1:30)
# ElbowPlot(pbmc,ndims = 30)
# dev.off()
# # find how many 15 cluster
# hms_neighbor<- FindNeighbors(pbmc, dims = 1:30)
hms_neighbor<- FindNeighbors(Duct_epithelial, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
pdf("./data/output/Duct_epithelial_cell_clustree_聚类树.pdf",width = 8,height = 6)
clustree(obj)
dev.off()


set.seed("123")
####resolution = 1.2####
hms_cluster <- FindClusters(hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
#UMAP可视化
pdf("./data/output/Duct_epithelial_cell_未标注_DimPlot.pdf",width = 7,height = 6)
DimPlot(hms_cluster, reduction = "umap",group.by = "tech",label = TRUE, pt.size = 0.5)
DimPlot(hms_cluster, reduction = "umap",split.by = "tech",label = TRUE, pt.size = 0.5)

highlight_clusters <- c("0","1", "2","3","4","5","6","7","8","9",
                        "10","11", "12","13","14","15","16","17","18","19",
                        "20") # 替换为您想要突出显示的亚群名称
highlight_colors <- c("white","white","white","white","white","white","white","white","white", "blue", 
                      "white","white","white","white","white","white","white","white","white","white",
                      "white") # 替换为对应的颜色
DimPlot(hms_cluster, reduction = "umap", split.by = "tech",label = TRUE, cols = highlight_clusters, cols.highlight  = highlight_colors)

dev.off()

#观察比例
table(hms_cluster@meta.data$seurat_clusters)
Cellratio <- prop.table(table(Idents(hms_cluster),hms_cluster@meta.data$tech) , margin = 2)#计算各组样本不同细胞群比例
Cellratio
# Normal 0 1 2 4 
# Normal       Tumor
# 0  0.210046512 0.019350074
# 1  0.168372093 0.020531758
# 2  0.153302326 0.025997046
# 3  0.018418605 0.130723781
# 4  0.155906977 0.020236337
# 5  0.065302326 0.083013294
# 6  0.018232558 0.090841950
# 7  0.047255814 0.067060561
# 8  0.020837209 0.076809453
# 9  0.059348837 0.044608567
# 10 0.026418605 0.063515510
# 11 0.003720930 0.080206795
# 12 0.006697674 0.077104874
# 13 0.011720930 0.062186115
# 14 0.026790698 0.029985229
# 15 0.005767442 0.041506647
# 16 0.000744186 0.023929099
# 17 0.000000000 0.018611521
# 18 0.000744186 0.013146233
# 19 0.000000000 0.005465288
# 20 0.000372093 0.005169867
# save
saveRDS(hms_cluster, file = "./data/temp/Duct_epithelial_cell_test_1.0.rds")

# cellmarkers用于后续标注
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/temp/Duct_epithelial_cell_cellMarkers.txt",sep="\t",row.names=F,quote=F)
# 挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)
write.csv(file="./data/temp/Duct_epithelial_cell_cellmarkers_top20.csv",top20)
# 整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/Duct_epithelial_cell_cellmarkers_top20_genes_1.2.csv",top20_table,row.names=F)


###2.2 亚群注释####
hms_cluster<-readRDS("./data/temp2/Duct_epithelial_cell_test_1.0.rds")
DimPlot(hms_cluster, reduction = "umap", label = TRUE)

#细胞及细胞中基因与RNA数量
slotNames(hms_cluster)
hms_cluster@assays
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)

#1.直接注释（未采用）⚠️
# new.cluster.ids <- c("group1", 	"group2", 	"group1", 	"group3", 	"group2", #4
#                      "group4", 	"group8", 	"group3", 	"group9", 	"group5", #9
#                      "group9", 	"group8", 	"group6", 	"group6", 	"group7", #14	
#                      "group7", 	"group3", 	"group9", 	"group9", 	"group3", #19
#                      "group8")
# names(new.cluster.ids) <- levels(hms_cluster)
# hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
# DimPlot(hms_cluster_id, reduction = "umap", split.by="tech",label = TRUE)

#2.注释前去掉17-20这些小亚群，美观⭐
hms_cluster1 <- subset(hms_cluster,idents=c(17:20),invert=TRUE)
DimPlot(hms_cluster1, reduction = "umap", label = TRUE)
new.cluster.ids <- c("Group_1", 	"Group_2", 	"Group_1", 	"Group_4", 	"Group_2", #4
                     "Group_3", 	"Group_5", 	"Group_4", 	"Group_7", 	"Group_6", #9
                     "Group_4", 	"Group_6", 	"Group_9", 	"Group_9", 	"Group_8", #14	
                     "Group_8", 	"Group_4")
names(new.cluster.ids) <- levels(hms_cluster1)
hms_cluster_id<- RenameIdents(hms_cluster1, new.cluster.ids)
#All
DimPlot(hms_cluster_id, reduction = "umap", label = TRUE)
# Normal_Tumor
DimPlot(hms_cluster_id, reduction = "umap", split.by="tech",label = TRUE)

# levels()修改ident的顺序，解决后期作图排序
levels(hms_cluster_id)
levels(hms_cluster_id) <- c("Group_1", 	"Group_2", 	"Group_3", 	"Group_4","Group_5", "Group_6", 	"Group_7", 	"Group_8", 	"Group_9")
levels(hms_cluster_id)

saveRDS(hms_cluster_id, file = "./data/temp/Duct_epithelial_cell_cluster_id_test.rds")


#整合后的cellmarker 20240918补充
hms_cluster <- readRDS("./data/temp/Duct_epithelial_cell_cluster_id_test.rds")
# cellmarkers用于后续标注
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/temp/Duct_epithelial_cell_cellMarkers_20240918补.txt",sep="\t",row.names=F,quote=F)
# 挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)
write.csv(file="./data/temp/Duct_epithelial_cell_cellmarkers_top20_20240918补.csv",top20)
# 整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/Duct_epithelial_cell_cellmarkers_top20_genes_1.2_20240918补.csv",top20_table,row.names=F)





###——UMAP出图####
#All
p1 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE)
p1
# Normal_Tumor
p2 <- DimPlot(hms_cluster_id, reduction = "umap", split.by="tech",label = TRUE)
p2
pdf("./data/output/Duct_epithelial_cell_sub_分群标注ALL.pdf",width =8 ,height = 6)
p1
dev.off()
pdf("./data/output/Duct_epithelial_cell_sub_分群标注Normal_Tumor.pdf",width =19 ,height = 10)
p2
dev.off()



###2.3 注释验证####
#绘图准备： 读入+亚群提取
hms_cluster_id<-readRDS("./data/temp/Duct_epithelial_cell_cluster_id_test.rds")
levels(hms_cluster_id)


####1.3.1 癌基因_Dimplot⚠️未使用####
#单独的基因分布,可以选择不同基因显示，按照需要作图
#可能后面需要写个循环出图
p1 <- FeaturePlot(hms_cluster_id,reduction ="umap", split.by="tech",features =c("TOP2A"))
p1
p2 <- FeaturePlot(hms_cluster_id,reduction ="umap", split.by="tech",features =c("FXYD2"))
p2
pdf("./data/output/Duct_epithelial_cell_sub_癌基因_小提琴图.pdf",width =7 ,height = 3)
p1
p2
dev.off()
####1.3.2 癌基因_VlnPlot####
#癌基因
p1 <- VlnPlot(hms_cluster_id, features = c("FXYD2"),pt.size = 0)
p1
p2 <- VlnPlot(hms_cluster_id, features = c("AMBP"),pt.size = 0)
p2

p3 <- VlnPlot(hms_cluster_id, features = c("REG3A"),pt.size = 0)
p3
p4 <- VlnPlot(hms_cluster_id, features = c("TOP2A"),pt.size = 0)
p4

pdf("./data/output/Duct_epithelial_cell_sub_癌基因_小提琴图.pdf",width =7 ,height = 3)
p1|p2
p3|p4
dev.off()


####1.3.3 DotPlot: 气泡热图####
#⚠️T亚群注释的并不清晰
#⚠️⚠️⚠️ene还要重新选的,可以再加gene
markers1 <- c(
  "FSTL5"	, "SBSPON"	, "FGA"	, 
  "KCNH6"	, "CARMN"	, "CFC1"	,
  "STRA6"	, "SHC4"	, "CSF2",
  "IGFL3"	, "FCER1A"	, "FAM92B"	, 
  "SCEL"	, "SERPINB2"	, "MUC16"	, 
  "CSPG4" ,"SLC14A1"		,"MIA",
  "SPINK1"	, "PIGR"	, "CLDN18"	,
  "TOP2A"	, "ARHGAP11A"	, "CEP55"	, 
  "MYBPC1"	, "PGC"	, "ANKRD36C")
  
p7 <-DotPlot(hms_cluster_id, #"apCAF",
             features = markers1, 
             dot.scale = 6) + RotatedAxis()
p7

pdf("./data/output/DEC_marker_气泡热图.pdf.pdf",width = 13,height = 8)
p7
dev.off()

####1.3.4 RidgePlot：Cellmarker表达量-Normal Vs Tumo ####
#Cellmarker表达量:对比Normal Vs Tumo
p9 <- RidgePlot(sub_cluster_id, 
                features = c(
                  "HLA-DRA","CD74","HLA-DRB1",# apCAF
                  "FBLN1", "PLA2G2A","HSPA1A", # iCAF 
                  "FN1", "COL1A2", "COL1A1", # mCAF
                  "CST1", "CST2", "MYH11") , # myCAF 
                cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())
p9
pdf("./data/output/Duct_epithelial_cell_sub_标注验证6_cellmarkers表达量对比.pdf",width = 10,height = 10)
p9
dev.off()

####1.3.5 DoHeatmap：全亚群-热图####
markers.to.plot <-c(
  "HLA-DRA","CD74","HLA-DRB1",# apCAF
  "FBLN1", "PLA2G2A","HSPA1A", # iCAF 
  "FN1", "COL1A2", "COL1A1", # mCAF
  "CST1", "CST2", "IGFBP2") # myCAF 
p10 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
                 features = markers.to.plot,
                 size=5,
                 group.bar.height=0.03)+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p10
T_sub_heatmap<-p10$data
write.table(T_sub_heatmap,"./data/temp/Duct_epithelial_cell_sub_heatmap.csv")

p11 <- DoHeatmap(subset(sub_cluster_id,downsample=50000),
                 features = markers.to.plot,
                 size=5,
                 group.by="celltype")+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p11
T_sub_celltype_heatmap<-p11$data
write.table(T_sub_celltype_heatmap,"./data/temp/Duct_epithelial_cell_sub_celltype_heatmap.csv")

pdf("./data/output/Duct_epithelial_cell_sub_标注验证7_Heatmap_genes_celltype.pdf",width = 16,height = 13)
p10
p11
dev.off()







